function [z m s] = gaussprod(m1, s1, m2, s2)

if nargin < 3
    m2 = m1;
    s2 = s1;
end

s = inv(inv(s1) + inv(s2));
m = s*(inv(s1)*m1' + inv(s2)*m2');
m = m';

z = 1/(sqrt(det(2*pi*(s1+s2))));
z = z * exp(-0.5*(m1-m2)*inv(s1+s2)*(m1-m2)');
